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The Onsager rule determines the frequencies of quantum oscillations in magnetic fields. We 
show that this rule remains intact to an excellent approximation in the mixed-vortex state of the 
underdoped cuprates even though the Landau level index n may be fairly low, n ~ 10. The models 
we consider are fairly general, consisting of a variety of density wave states combined with d-wave 
superconductivity within a mean held theory. Vortices are introduced as quenched disorder and 
averaged over many realizations, which can be considered as snapshots of a vortex liquid state. We 
also show that the oscillations ride on top of a held independent density of states, p{B), for higher 
helds. This feature appears to be consistent with recent specihc heat measurements [C. Marcenat, 
et al. Nature Comm. 6, 7927 (2015)]. At lower helds we model the system as an ordered vortex 
lattice, and show that its density of states follows a dependence p{B) oc %/B in agreement with the 
semiclassical results [G. E. Volovik, JETP Lett. 58, 469 (1993)]. 


I. INTRODUCTION 

A breakthrough in the area of cuprate supercon¬ 
ductivity is the observation of quantum oscillations in 
cuprates m m- In these experiments a strong mag¬ 
netic field is applied to suppress the superconductivity, 
which most likely reveals the ground state [5] without 
superconductivity. However, the understanding of this 
“normal state” may be a crucial ingredient in the the¬ 
ory high temperature superconductivity. Standing in the 
way are at least two important issues: (1) Does the quan¬ 
tum oscillation frequencies substantially deviate from the 
classic Onsager rule for which the oscillation frequency 
F = {hc/2TTe)A{ep), where ^(cf) is equal to the ex¬ 
tremal Fermi surface area normal to the magnetic field? 
If so, it would lead to considerable uncertainty in the in¬ 
terpretation of the experiments. (2) Do the oscillations 
ride on top of a magnetic field dependence of the density 
of states(DOS) p{B) ^ y/Bl [1] If so, it might indicate 
the presence of superconducting fluctuations even in high 
magnetic fields at zero temperature, T = 0, from an ex¬ 
trapolation of a result of Volovik, which is supposed 
to be asymptotically true as H —>■ 0. Therefore high field 
behavior requires careful analyses. Quantum oscillations 
require the existence of Landau levels. If this is true, 
they might indicate the existence of normal Fermi liquid 
quasiparticles [H]. 

It has been argued from a theoretical analysis that the 
Onsager rule could be violated by as much as 30% [7] . We 
find that under reasonable set of parameters, to be de¬ 
fined below, the violation is miniscule, ~ 10“"^. Even for 
extreme situations discussed in Ref. [TJ it is less than 2%. 
If we are correct, one can use the Onsager rule to inter¬ 
pret the experiments with impunity. The second encour¬ 
aging result is that p{B) saturates in the regime where 
oscillations are present. We interpret this to mean that 
there are generically no superconducting fluctuations in 
high fields. A recent specific heat measurement [8] shows 
that the specific heat indeed saturates at high fields, sig¬ 


nifying that the normal state is achieved. 

To put our paper in the context, note that in conven¬ 
tional s-wave superconductors, previous work has shown 
that for higher Landau level indices, and within coher¬ 
ent potential approximation, vortices mainly damp the 
oscillation amplitude, but the shift in the oscillation 
frequency [9l llOj is negligible; however, for d-wave un¬ 
derdoped cuprate superconductors with small coherence 
length and high fields with Landau level indices ~ 10 
this calculation should not hold [7]. A more recent semi¬ 
classical analysis based on an ansatz of gaussian phase 
fluctuations of the d—wave pairing m indicates that 
the oscillation frequency is unchanged, as here. However, 
relatively undamped quantum oscillations riding on top 
of '/H was found in this dynamic Gaussian ansatz that 
does not account for vortices, which must necessarily be 
present, as in Ref. [7], and the branch cuts introduced by 
the vortices must also be taken into account. 

We consider the vortices explicitly in the Bogoliubov- 
de Gennes (BdG) Hamiltonian, as in Ref. [7], and model 
the vortex liquid state as quenched, randomly distributed 
vortices, paying special attention to branch cuts. There 
are other important differences as well, as we shall dis¬ 
cuss below. To have a complete picture, we also consider 
the low field regime where the quantum oscillations dis¬ 
appear. In this regime the vortices arrange themselves 
into a vortex solid state and should be modeled as an 
ordered lattice instead. We compute the DOS of such 
vortex lattices explicitly and find that p{B) oc '/B in the 
asymptotically low field limit, consistent with Volovik’s 
semiclassical analysis [5]. 

In Section H we define the model Hamiltonian that in¬ 
cludes d-wave superconducting order parameter as well as 
a variety of density wave states. Our numerical method, 
the recursive Green function method adapted for the 
present problem is discussed in Sec. HI. The results are 
discussed in Sec. IV and Sec. V contains discussion. 
There are three appendices. 
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II. THE MODEL HAMILTONIANS 

The starting point is the Bogoliubov-de Gennes (BdG) 
Hamiltonian 


n = 


H - fi Aij \ 

Al, + 


( 1 ) 


defined on a square lattice. Here is the chemical po¬ 
tential. H is the Hamiltonian that describes the normal 
state electrons; while the off-diagonal pairing term 
defines the superconducting order parameter. For sim¬ 
plicity, we ignore self consistency, as we believe that it 
cannot change the major striking conclusions. 


A. The diagonal component H 

Besides the hopping parameters, the normal state 
Hamiltonian H contains a variety of mean field order pa¬ 
rameters defined below. Although many different orders 
are suggested to explain the normal state of the high- 
Tc superconductivity, for our purposes it is sufficient to 
consider three different types: a period-2 d—density wave 
(DDW) [SI Uni [13], a bi-directional charge density wave 
(CDW), and a period—8 DDW model [T3|. We believe 
that our major conclusions in this paper do not depend 
on the nature of the density wave that is responsible for 
Fermi surface reconstruction. The DDW is argued to be 
able to account for many features of quantum oscillations, 
as well as the pseudogap state [Hds] in the cuprates. 
Among the many different versions of density waves of 
higher angular momentum nn, the simplest period-2 
singlet DDW, also the same as staggered flux state in 
Ref. |7], and also a period—8 DDW order, proposed by 
us previously to explain quantum oscillations [13], have 
been chosen here for illustration. 

Recently a bi-directional CDW has been observed 
ubiquitously in the underdoped cuprates [IMI]- It 
has ordering wavevectors Qi « ^(0.31,0) and Q 2 « 
— (0,0.31), which are incommensurate. This order has 
also been used to explain the Fermi surface reconstruc¬ 
tions and quantum oscillation experiments [221 [23], al¬ 
though a recent numerical work [24] has demonstrated 
that the strict incommensurability of the CDW can de¬ 
stroy strict quantum oscillations completely. For the pur¬ 
pose of illustration we chose, instead, commensurate vec¬ 
tors Qi = ^(|,0) and Q 2 = ^(0,1). 

Therefore without the magnetic field, B, H is given by 

H =-tJ2 4.Cr, +t' Y. cjTr, 

-t” Y. + h.c.-f Hd.w., (2) 

((dj)» 

where t,t', and t" are the 1st, 2nd and the 3rd nearest 
neighbor hopping parameters respectively. Hd.w. is var¬ 
ious density wave orders specified below. The external 
uniform magnetic field B = Bz is included into H via the 


Peierls substitution: cJ.Crj => exp[—J))'’A • dlJcJ.Crj 
with the vector potential A. = B xy chosen, for simplic¬ 
ity, in the Landau gauge. 

1. T-wo-fold DDW order 

where S = x,y denote the two nearest neighbors. 
r]s = 1 ioT S = X while = — 1 for <5 = y indicates 
the DDW order has a local d—wave symmetry. 

2. Bi-directional CDW order 

Hd.w. =Vc^r]s {cos[Qi • (r^ -F 5/2)] 

-h cos[Q 2 • (ri -h 5/2)]} cl__f_gCri. (4) 

Again 775 = ±1 is the local d—wave symmetry factor 
of the CDW order. 

3. Period-8 DDW order 

-^d.w. — 'y ^ Cj^Ckd_Q -t- V); Cj^Ck_|_ 2 Q -l- h.c. . (5) 

k 

where the iGk term is the period —8 DDW order 

^ <^kyk-i-Q ^Ck' ^k.k^d-Q (b) 

Here Gk = (Wk - B^k+Q)/2 with the DDW 
gap Wk = ^{coskx — cosky), and the ordering 
wavevector is Q = (|^, ^). The 14 term in Eq. ([^ 
represents a period—4 unidirectional CDW with an 
ordering wavevector 2Q, which is consistent with 
the symmetry of the period —8 DDW order. Notice 
that this CDW is different from the bi-directional 
CDW we considered in the previous section. Exper¬ 
imentally whether the observed CDW in cuprates 
is unidirectional or bi-directional is still not fully 
resolved. 

Fourier transformed to the real space, the Hamil¬ 
tonian Hd.w. becomes 


Hd. 


W. 


^,Wo . Q-(r-r') . Q-(r + r') 
2 2 2 


^ 5i-.x' — ax 5x.r'd-ay 

+2 14 ^ cos[2 Q • r] clcr- 



r 


In the above, Wq controls the overall magnitude of 
the period-8 DDW order parameter, i sin ^ ^ 

indicates that the order is a current, sin ^ 

shows that the magnitude of this order is modu¬ 
lated with a wavevector Q, and the last factor in 
the curly bracket {...} explicitly exhibits its local 
d—wave symmetry. If Q = (tt/o, tt/o), then the 
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order reduces to the familiar two-fold DDW order. 
In that case | sin ^ | = 1 and the order pa¬ 

rameter magnitude is a constant. The last term 
in Eq. 0 gives the 2Q charge modulation. Note 
that this CDW is defined on sites, differing from 
the bi-directional CDW defined on bonds. 


B. The off-diagonal component Aij 

The off-diagonal pairing term A.y in the BdG Hamil¬ 
tonian is defined on each bond connecting two nearest 
neighboring sites i and j. Aij = rjij, where 

Tjij = -|-1 if the bond is along x—direction and rjij = —1 
if it is along j/—direction so that A^^ has a local d—wave 
symmetry. The pairing amplitude is taken to be 


with m* = 2m and e* = — 2e are the mass and the charge 
of the Cooper pairs respectively. The path for this inte¬ 
gral is chosen such as to avoid the branch cuts of all the 
vortices so that the phase field (j)i is single valued on every 
site, as illustrated in Fig. 



|A,,| = A 


7’eff 


( 8 ) 


where A is the pairing amplitude far away from any vor¬ 
tex center. ^ is the vortex core size. In our calculation 
^ = 5a is adopted, where a is the lattice spacing. In the 
presence of a single vortex, reS in the above is simply 
the distance from the center of our bond to the 

center of that vortex. While in the presence of multiple 
vortices, following the ansatz used in Ref. [7] we choose 
(r^)^ “ where r„ is the distance from the bond 

center to the nth vortex center and g > 0 is some real 
number. 

In this ansatz, VeS is a monotonic increasing function of 
the parameter q for a given vortex configuration. There¬ 
fore if q is large, the calculated rgg as well as |Aij| is 
also larger, which means the vortex scattering is stronger. 
However our conclusions do not depend on the different 
choices of q (for more details see the appendix section]^. 
Therefore in this paper, if not specified otherwise, q = 2 
will be chosen. 

The bond phase variable Oij contains the information 
of our quenched random vortex configuration, but for 
the purpose of our calculation we need the site phase 
variables. We use the ansatz for 0^- given in Refs. P51[2S] 


e 


i9i. 


4>i 

= e* 2 sgn[cos 



(9) 


where 4>i is the pairing order parameter phase field de¬ 
fined on a site. In the above, without the “sgn[...]” factor 
Oij is simply the arithmetic mean of (j)i and (pj. However 
using Oij = is not enough because whenever the 

bond ij crosses a vortex branch cut, the phase factor 
will be incorrect and different from the correct one 
by a minus sign. This can be corrected by the additional 
“sgn[...]” factor(see the appendix section 0. 

Then pi can be further computed from the superfluid 
velocity field Vs(ri) by 

<>.-* = CI=^ + liA(r)].dl (10) 


FIG. 1. Illustrations of the vortices (circle) on the lattice 
(dashed lines). The arrows show the path of the integral we 
have chosen in defining our phase field 4>(r). To make this 
phase definite, the branch cuts of all the vortices are chosen to 
extend from the vortex center to the positive infinity [x = oo), 
represented by the magenta horizontal lines. 


We still need to compute the superfiuid velocity Vg (r). 
This can be done by following Ref. m 


mVg(r) = —mh J 


(27r)2 -h A-2 ^ 

^ ^ r>. 


( 11 ) 


where m is the electron mass, A is the penetration depth, 
and R„ gives the nth random vortex position. In this 
integrand, because k x i is odd in k, only the imaginary 
part of e*k (r-R„) survive after the integration, so the 
whole expression on the right hand side becomes real. We 
also make an approximation A = oo so that we can ignore 
the term in the denominator. This is equivalent to 
replacing the magnetic field B(r) by its spatial average, 
which is equal to the external magnetic field B = Bz. 
It is a good approximation when B 3> Hci, where Hd 
is the lower critical field. This condition is well satisfied 
in the quantum oscillation experiments of cuprates. Also 
this approximation is consistent with our initial choice of 
the vector potential A = Bxy, given completely by the 
applied external field B. 

For our square lattice calculation we discretize the 
above k integral and choose 27r/^ as its upper cutoff, 
since the vortex is only well defined over a length scale 
larger than the vortex core size Therefore in the limit 
A ^ ^ > a, Vs(r) can be rewritten as follows 

^ LMa^ ^ ^2^sm[k.(r-R„)]. 

(fcx,fc») 

( 12 ) 


In this summation -|- ...., ^ - ll, 

fcy = - ^, - ^ -b ^ ^^. The prime super¬ 

script in the summation means the point (fc^,, ky) = (0, 0) 
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is excluded to be consistent with our approximation 
A = oo. 


III. THE RECURSIVE GREEN FUNCTION 


Given the BdG Hamiltonian % defined above, we use 
the recursive Green’s function method [21] to compute 
the local DOS(LDOS). We attach our central system, 
which has a lattice size L x M, to two semi-infinite 
leads in the ix directions. The leads are normal met¬ 
als described by t,t',t'' only. Then we can compute 
the retarded Green’s function Gi{j,j'-,E -|- i6) at an en¬ 
ergy E for the ith principal layer(see the appendix sec¬ 
tion E- Here each ith principal layer contains two adja¬ 
cent columns of the original square lattice sites. So there 
are L/2 principal layers and each of them contains 2M 
number of sites. Therefore E + iS) is a 4M x AM 

matrix, with j,j' = 1,2, ...,AM, because it has both an 
electron part and a hole part. In calculating the LDOS 
at the jth site of the Hh layer only the imaginary part 
of the jth diagonal element in the electron part of Gi 
is included. This is equivalent to treating the random 
vortices as some off-diagonal scattering centers for the 
normal state electrons. To see smooth oscillations of the 
DOS we also average the calculated LDOS over different 
sites and realizations of uncorrelated vortices. In other 
words the quantity of our central interest is 


p{B) = 



LI2 2M ^ 

^ ^-)Im Gi(j, j; 0 -b z d) 

i=i j=i 


(13) 


where the angular brackets denote average over indepen¬ 
dent vortex realizations. In the Green’s function we have 
already set the energy to the chemical potential E = 0. 
For all the numerical results presented in the following, 
an infinitesimal energy broadening S = 0.005t will be cho¬ 
sen, if not specified otherwise, and the periodic boundary 
condition is imposed in the ?/—direction. 


IV. RESULTS 

A. The Onsager rule for quantum oscillation 
frequencies 

1. The two-fold DDW order case 

With the parameters: t = l,f' = 0.30 = 

t'/9.0,p. = —0.88071, Wo = 0.26 f, 14 = 0, the hole dop¬ 
ing level is p « 11%. Without vortices we can diagonalize 
the Hamilontian E[ in the momentum space and obtain 
the normal state Fermi surface. This Fermi surface con¬ 
sists of two closed orbits, see the inset of Fig. The 
bigger one centered around the node point (f, f) is hole 
like. It has an area ~ 3.47%. This corresponds 

to an oscillation frequency Fh = = 966T 


from the Onsager relation, where = ft,c/2e is the 
fundamental flux quanta and the two lattice spacings 
= 3.82A X 3.89A are chosen for YBGO. At the antin- 
odal point (0, tt) there is an electron pocket with an area 
( 2 ^fay ~ 1-9%, corresponding to a frequency E^ = 525T 
(electron). We should notice that the fast oscillation Fh 
(hole) is not observed in the experiments in cuprates. 
This problem can be resolved if we consider a period—8 
DDW model |T3|; see below. 

We compute the p{B) as a function of the inverse of 
the magnetic field 1/F in the presence of various A. In 
these calculations, the number of vortices are chosen such 
that the total magnetic flux is equal to $ = BLMa^. 
From the oscillatory part of p{B) we perform Fast Fourier 
Transform(FFT) to get the spectrum. The result is 
shown in Fig. |2d[ In this spectrum the two oscillation 
frequency Fe = 525T and Fh = 966T calculated from the 
normal state Fermi surface areas via the Onsager relation 
are also shown by the two vertical dashed lines. We see 
clearly that as we increase A the oscillation amplitudes 
are damped. However, remarkably, the oscillation fre¬ 
quencies remain the same within numerical errors. Thus, 
even in the presence of vortices, the Onsager rule still 
holds to an excellent approximation. 


2. The bi-directional CDW order case 


We choose the following parameters: t = l,t' = 
0.2t, t" = t'/8, Vc = 0.12t, p, = —0.73t so that we can pro¬ 
duce the right oscillation frequencies that are observed 
in experiments. The hole doping level is p « 11%. The 
Fermi surface of the normal state is plotted in the inset 
of Fig. 2b (open orbits are not shown for clarity). There 
are two closed Fermi surface sheets. Gentered around 
the point (f, f) and other symmetry related positions 
there are diamond shaped electron pockets, highlighted 
in orange. This pocket has an area ( 2 ^ 0 )^ ~ 1.9%. It 
corresponds to a frequency Fe = 529T from the Onsager 
relation. Besides this electron pocket, there is an oval 
shaped hole pocket centered around (|, ^), highlighted 
in blue. The area of this hole pocket is = 0.33%. 

This corresponds to an oscillation frequency Fh = 92T. 


The oscillation spectrum of the piB) is shown in 
Fig. From the spectrum we see that when the vortex 
scattering is absent, A = 0, the oscillation amplitudes 
peak at the two frequencies Fe,Fh, as denoted by the 
two vertical dashed lines. These results agree with our 
Fermi surface calculation, as we expected. When the 
vortices are included the oscillation amplitude is gradu¬ 
ally damped as the vortex scattering strength is increased 
by increasing A. However whenever the oscillation fre¬ 
quency can be clearly resolved, we see that their positions 
do not change with A. Again this means that the On¬ 
sager rule survives in the presence of vortex scattering. 








5 






✓ V 


\ ✓ 




Ea 


- 3 - 2-10 1 2 3 



D ^ a 


^ ^ 

ay Sr Vsk 


Q Q 


0.0 0.5 1.0 1.5 2.0 2.5 3.0 


(a) FS for the two-fold DDW 



(d) FFT spectrum for the two-fold DDW 


(b) FS for the bi-directional CDW 


(c) FS for the period-8 DDW 




(e) FFT spectrum for the bi-directional (f) FFT spectrum for the period-8 DDW 
CDW 


FIG. 2. The upper panel shows the plots of Fermi surfaces(FS), with electron pockets shaded in orange and hole pockets in blue. 
In Fig. |2a| the area enclosed by the dashed lines gives the reduced Brillouin zone. In Fig. |2b| the open orbits are not shown for 
clarity. And the dashed lines denote the positions kya = tt/S, 27r/3. In the lower panel we show the corresponding oscillation 
FFT spectrums for various values of A. The two vertical dashed lines in each plot denote the two fundamental oscillation 
frequencies calculated from the Fermi surface area via the Onsager relation. They are Fe = 525T, Fh = 966T in Fig. |2d| 
Fe = 529T, Fh = 92T in Fig.[^ and Fe = 523T, Fh = 159T in Fig.[^ Other parameters used are L — 1000, M = 100, S = O.OOSt 
in Fig. L = 1000, M = 102, 5 = 0.005t in Fig. and L = 1000, M = 200 ,5 = 0.002t in Fig. 


3. The periods DDW order case 

In this subsection we present our quantum oscillation 
results for the period—8 DDW model. In this model the 
period—8 stripe DDW order is considered as the major 
driving force behind the Fermi surface reconstructions; 
while a much weaker unidirectional period—4 CDW is 
included as a subsidiary order. 

We choose the parameter set t' = 0M,t” = 

t'l2.0,Wo = 0.70t,Vc = O.OSt,^ = —0.70t and estimate 
the hole doping level to be p « 11.2%. We also obtain 
a Fermi surface similar to the one we had in Ref.|18j. 
It has a large pocket of electron like with a frequency 
Fe = 523T, a smaller pocket of hole like with a frequency 
Fh = 159T, and also some open orbits which do not con¬ 
tribute to quantum oscillations. 

The corresponding oscillation spectrum is presented in 
Fig. 1^ where we see the oscillation amplitude decreases 
as we increase A, however, the frequencies do not change 
with A. In other words the presence of vortex scattering 


does not alter the oscillation frequencies. 

The observations here, combined with the other two 
cases, strongly suggest that the Onsager’s relation being 
intact in the presence of vortex scattering is generic and 
independent of the order parameters that reconstruct the 
Fermi surface. 


B. The Density of states at high fields 


In the above we have examined the effects of random 
vortex scattering on the quantum oscillations. Now we 
give an overview of the B dependence of the DOS for 
fields B > lOT, at a representative value of A = O.lt. 
At lower fields, the vortex liquid model is not valid any 
more, since vortices should order into a solid instead. 
Therefore we should use a vortex lattice to model such a 
state. In the following we focus on the high field regime 
first and defer our vortex lattice discussions for the low 


field regime to the later section IV C 
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1. The period-2 DDW order case 


In Fig. 3a we plot p(i3)/p„(0) as a function of the field 
B for the two-fold DDW order case, where Pn(0) is the 
normal state DOS at zero field. In the following, the nor¬ 
mal state should be understood as a state, which does not 
have any superconductivity but can have a particle-hole 
density wave order. And all the DOS value calculated is 
for one electron in a single CuO plane, without includ¬ 
ing the spin degeneracy. From Fig. [3^ we see that as 
B decreases, the DOS oscillation gets suppressed gradu¬ 
ally. This is because the orbital quantization of electrons 
becomes dominated by the vortex scattering. 

A noticeable feature of this plot is that when the field 
becomes large, the oscillation of p{B) in \/B gradually 
develops on top of a constant background. This constant 
background value of p{B) is suppressed from the normal 
state DOS Pn(0). The size of this suppression depends 
on the vortex scattering strength. For the parameters 
used in Fig. [^it is ~ 15%. This constant background 
of p{B) is different from the previous results obtained in 
Fig. 3(b) of the Ref. [TT] in the absence of vortices. 


2. The bi-directional CDW order case 


The constant background of the DOS oscillation is not 
restricted to the two-fold DDW order case. As we can 


see in Fig. 3b for the bi-directional CDW order, the p{B) 


oscillation background is again a constant at high fields. 


3. The period—% DDW order case 

We also confirm this constant p{B) background feature 
in the oscillation regime for the period—8 DDW order 
case in Fig. 

Therefore we can conclude that the high field p{B) 
oscillation background being a constant is generic. 


C. Vortex solid at low fields 

Now we move on to the low field regime. In this regime 
when the field is low enough, the vortices order into a lat¬ 
tice. Whether the lattice is square or triangular requires 
a self-consistent computation of the system’s free energy, 
which is far beyond the scope of this paper. Instead we 
simply take a square lattice for illustration. But none of 
the following qualitative features should depend on the 
vortex lattice type. 


preserved along the transverse direction, we require the 
vortex lattice to be commensurate with our original CuO 
lattice, as schematically shown in Fig. Namely, if the 
vortex lattice spacings are [lx, ly), and the corresponding 
vortex lattice size is {Nx, Ny), we require that the original 
CuO lattice size {L,M) satisfies L = Nxlx,M = Nyly. 
For a particular value of {L,M), this restricts the pos¬ 
sible values of {lx,ly) and also the possible values of 
the magnetic field, because the vortex lattice spacings 
{lx,ly) are connected to the magnetic flux density via 
B — where d)g = /ic/2e is the fundamental 

flux quanta. In our following calculation we pick a partic¬ 
ular value of the system size (L, M), find all the possible 
compatible values of the vortex lattice spacings lx = ^y, 
and then for each of them calculate the magnetic field B 
as well as the corresponding DOS. 

However, we should calculate the DOS of the Bogoli- 
ubov quasiparticles instead of the electrons, because the 
system is far from being in a normal state in such a low 
field regime. Therefore now p{B) is computed from the 
following formula instead 

L/2 am 

P{B) = Gi{j, j-,{) + i 5). (14) 

i=l 3 = 1 


The major differences here from the one we used in our 
quantum oscillation calculations are: (1) the summation 
of the Green’s function’s diagonal matrix elements in¬ 
cludes both the electron part and the hole part: j runs 
from j = 1 to j = AM instead of j = 2M; (2) there is no 
averaging over different vortices configurations because 
the vortex lattice is ordered; (3) an additional prefactor 
of 1/2 is added to avoid double counting of degrees of 
freedoms. 

For such a vortex lattice calculation, the summation 
over different vortex positions in the superfluid velocity 
calculation in Eq. (121 can be done exactly by using 


E' 


,ik(r 


= NxNy 


E ' 

("l,n2) 


(15) 


where = (^y/y, x?) is a reciprocal Bragg vector 

of the square vortex lattice, with ni,n 2 € Z. Then the 
Eq. (12) of Vs becomes 


mv, = Tih 


1 


Ixhi 


E ' Gni n2 ^ ^ 

^ sm[G, 


(ni.n2) 


\Gnuu,\^ 


(16) 


The summations of (ni, 712 ) are restricted to those values 
that satisfy 0 < < —, 0 < < —. Again 

the prime superscript in the summation means the point 
( 711 , 712 ) = (0)0) is excluded. 


1. Implementation of the square vortex lattice 

2. DOS numerical results 

To put the square vortex lattice onto our original CuO 

lattice so that each vortex sits at the CuO lattice pla- According to Volovik [3, for a -y 2 —wave vortex, the 
quette center and the periodic boundary condition is still major contribution to the low energy DOS comes from 
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(b) The bi-directional CDW order 


(c) The period-8 DDW order 


FIG. 3. The DOS p(B),normalized to the normal state DOS Pn(0) at zero field B — 0 for different cases. The estimated values 
of the normal state DOS are Pn(0) « 0.23states/t for the two-fold DDW, Pn(0) « 0.25 states/f for the bi-directional CDW, 
and Pn(0) « 0.18 states/f for the period-8 DDW case. The data is averaged over 108, 120,40 different vortices configuration 
realizations respectively. 
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"I" where Aa ~ 90T, if ^ = 5a and a « 3.83A are used. 
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^ I. F irst we consider a square vortex lattice without 

1 any other additional density wave order. We choose 

the band structure parameters to be t'/t = 0.3, t" = 


FIG. 4. Schematic diagram of a square vortex lattice. The 
dashed lines represent the original CuO lattice; while the full 
lines stand for the vortex lattice, with each vortex, repre¬ 
sented by the grey disks, sitting at the CuO plaquette center. 
And {lx,ly) are the vortex lattice spacings, in units of the 
original CuO square lattice spacing a. 


the extended states along the nodal direction. In his 
semiclassical analysis this contribution is computed from 
the Doppler shift of the quasiparticle energy. The conclu¬ 
sion is that the DOS for a single vortex is p{B) oc Ij'/B. 
In the limit that the number of vortices is proportional 
to S, which is not valid if B is near the lower critical 
field 77ci, multiplying it by the number of vortices gives 
p{B) oc c/B. Extrapolating this result to the high field 
regime and using the fact that near the upper critical 
field Hc 2 , p{B) should roughly recover the normal state 
DOS p„(0), he concluded that p{B)/pn{0) = Ky/BjH^, 
with K some constant of order unity. This type of analy¬ 
sis is applicable only in the small field limit in the sense 
that B ^ Hc2 so that each vortex is far apart from any 
others. This is exactly the field regime where the vortex 
solid state develops. In the following we compute the 
DOS for a d—wave vortex lattice, for the cases both with 
and without an additional particle-hole density wave or¬ 
der, and test them against Volovik’s results. For our 
following comparisons we slightly rewrite the above field 


t'/9.0,p, = —l.Olt so that the estimated normal 
state hole doping level p « 15% is at the optimal 
doping. The computed DOS is shown in Fig. 5a 
At low enough fields, all the data points follow 
the p{B)/pn{0) = 0.3c/B line, although there is 
some small scatter in the data, which comes from 
the finite size effects of our vortex lattice. This 
O.Sc/B corresponds to k « 3 in Eq. 0- To make 
a comparison with the specific heat measurements 
on YBC0123 at the optimal doping [2S]) we es¬ 
timate the field dependent electronic specific heat 
7 (i?) from our DOS p{B) as follows 




'In 


Pn(0) 


(18) 


The normal state specific heat can be estimated as 
7 „ = 4 ^k^pniO). Here the additional prefactor of 
4 comes from the spin degeneracy and the fact that 
one unit cell of YBC0123 contains two CuO planes. 
If we take t = 0.15eV, then 7 „ « 15.7mJ/mol • 
and 7 (i?) = Ac/B with the coefficient A « 
4.7mJ/mol • • T^/^. Compared with the exper¬ 
imental value of A « 0.9 mJ/mol • • T^/^ from 

Ref. [2^ , our numerical value is greater by a factor 
of about 5. This quantitative discrepancy is not 
significant given our approximations. In fact, it is 
quite reasonably consistent. 


























(a) The pure d-wave vortex lattice case without any 
other density wave order at the optimal hole doping 



(b) The coexistence of a vortex lattice and a two fold 
DDW order in the underdoped regime 


FIG. 5. The DOS of vortex solids for A = O.lf. In Fig. 5a Pn(0) « 0.25 states/t and in Fig. 5b Pn(0) ~ 0.23 states/t. 


2. Next we consider the coexistence of a square vor¬ 
tex lattice and an additional two-fold DDW order 
in the underdoped regime. The parameters are 
the same as those in our high field quantum os¬ 
cillation calculations: t'/t = 0.3, t" = t'/9.0,fi — 
—0.8807t, Wq = 0.26t, so the estimated normal 
state, with the DDW order but no superconductiv¬ 
ity, hole doping level is p « 11%. Fig. ^ shows the 
corresponding DOS results. The small field data 
follows p{B)/pn{0) = OA'/B, corresponding to a 
value of K « 4 in Eq. 


The above two values of k are consistent with the fact 
that in Volovik’s formula k is of order unity. Of course its 
precise value depends on the vortex lattice structure, on 
the slope of the gap near the gap node (in the current case 
both the parameters A and q), and also on the normal 
state band structure. 

From the above two scenarios we can conclude that 
irrespective of the existence of an additional density wave 
order, the DOS of a clean vortex lattice always scales as 
p{B) oc '/B in the low field limit. 


sity of states follows p{B) cx '/B in the asymptotically 
low field limit, in agreement with Volovik’s semiclassical 
predictions. However in contrast to the previous sugges¬ 
tion our results show that this small field limit does not 
extend to the high field oscillatory regime of the vortex 
liquid state. Instead when the oscillations can be re¬ 
solved, the non-oscillatory background of p{B) flattens 
out, and becomes field independent consistent with the 
more recent specific heat measurements [8j. 
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V. CONCLUSION 

In summary we have shown that in the quenched vor¬ 
tex liquid state the quantum oscillations in cuprates can 
survive at large magnetic fields. Although the oscillation 
amplitude can be heavily damped if the vortex scatter¬ 
ing is strong, the oscillation frequency is given by the 
Onsager rule to an excellent approximation. Of course, 
when the field is small the quantum oscillations are de¬ 
stroyed by the vortices and p{B) gets heavily suppressed 
due to the formation of Bogoliubov quasiparticles. When 
the field is small enough, a vortex solid state forms in¬ 
stead and it can be modeled by an ordered vortex lattice. 
We show the field dependence of the vortex lattice’s den¬ 


Appendix A: Bond phase field 6ij of Aij 

The phase field Oij is defined on the bond ij, which 
connects two nearest neighboring sites i and j. Therefore 
it is natural to use the phase fields (j)i and (j)j > on the site 
i and site j respectively, to define ^ However 

this definition does not guarantee that whenever a closed 
path encloses a vortex, Oij along that path will pick up 
a 27r phase as the vortex is winded once. Therefore this 
Oij can not give the correct vortices configuration. It is 
incorrect whenever a vortex branch cut is crossed. To 
see this clearly, we map the phase field 4 >i along a closed 
path that encloses a vortex onto a unit circle since ipi is 
defined only modulo 27r, as schematically shown in Fig.|^ 
In this figure, the blue arc segment corresponds to the 








9 


bond ij on the closed path. Therefore an appropriate 9ij 
should be equal to some value of the phase field on this 
segment. When the bond ij does not cross any branch 
cut, 9ij = is indeed on the blue segment and can be 

a good definition of 9ij, as illustrated in Fig. 6 a; however, 
if the bond ij crosses a branch cut, we see that , 
indicated by the red arrow in Fig. | 6 b[ is not on the blue 
segment and can not be an appropriate definition of 9ij. 


In this latter case, 9ij = _ 7 ^ instead can be a 

good definition, since it falls onto the blue arc segment, 
as indicated by the blue arrow in Fig. | 6 b| 




(a) Bond ij does not cross the (b) Bond ij crosses the branch 
branch cut: |0i — </>j| < f cut: \(f>i — (pjl > it 

FIG. 6 . The black dot in the center represents a vortex. The 
dashed line, extended to the infinity, is its branch cut. The 
blue arc segment corresponds to the bond ij on a closed path. 
The phases (pi, cpj are measured counter -clo ck wisely from the 


upper side of the branch cut. In Fig. 6 a the bond ij does 


not cross the branch cut, and {(pi + (pj)j2 is a good definition 
for 9ij ; however, if the bond ij crosses a branch cut, as in 
Fig. | 6 b[ then {(pi -\- <pj)l2 can not be a correct definition of 
Oij. Instead {(pi + (pj)/2 — -n gives an appropriate definition of 
da. 


Based on these two scenarios, a good definition of 9ij 
will be 


e*®- = sgn[cos (Al) 

This definition of 9ij guarantees that whenever the phase 
field (pi along a closed path crosses a branch cut once, 
the defined 9ij crosses the same branch cut once as well. 
When there are multiple vortices enclosed, we only need 
to linearly superpose the contributions from each vortex 
together to the field (pi and 9ij respectively. It is not 
difficult to see that the above definition of 9ij is still good 
in these cases. For our numerical calculation convenience, 
we rewrite the above definition of 9ij in a slightly different 
way 


e 


iO ij 


gi<Pi _|_ gi'Pj 
\(,ij>i _|_ gi(t>j\ 


(A2) 


Appendix B: Recursive Green’s function method 

The recursive Green’s function method studies a quasi- 
one-dimensional system, which is a square lattice with a 
very long axis of length La along the a;—direction and a 
shorter axis of width Ma along the y—direction in our 



FIG. 7. Schematic diagram of the recursive Green’s function 
calculation. The sites enclosed by the blue dashed rectan¬ 
gle define the ith principal layer. The exact Green’s func¬ 
tion Gi has two self-energy contributions from both the left 
semi-infinite stripe and and the right one. The left stripe is 
characterized by its surface Green’s function Gf_i with the 
{i — l)th layer its surface; while the right one is characterized 
by another surface Green’s function Gpj^i with the {i + l)th 
layer its surface. 


problem. The system can be built up recursively in the 
a:—direction by connecting many one-dimensional stripes 
together. Each stripe has a direct coupling only to its 
nearest neighboring ones. This property is essential for 
the recursion. In our Hamiltonian "H the third nearest 
neighbor hopping t" provides the farthest direct coupling 
along the x—direction. It connects two sites 2a apart. 
Therefore each stripe necessarily contains two columns 
of the square lattice sites so that the direct coupling ex¬ 
ists only between two adjacent stripes. The blue dashed 
rectangle in Fig. shows one such stripe. We define each 
of such stripes as a principal layer, so each layer contains 
2M sites. 

Our goal is to compute the diagonal matrix elements 
of the exact Green’s function G in order to get the DOS. 
For this purpose we first calculate Gi =< i\G\i > for 
each layer i. The ket \i > represents a state where the 
Bogoliubov quasiparticles are found in the fth principal 
layer. It has AM components, of which the first 2M ones 
give the electron part wavefunction, while the rest 2M 
ones define the hole part. Therefore Gi = Gi{j,j') is a 
AM X AM matrix with j,j' = 1,2,...,4M. For brevity 
we will suppress the matrix element indices hereafter, if 
there is no confusion. 

The exact Green’s function Gi can be computed by 
(for derivations see Ref. [55]) 


1-1 


^ _ r/^U ^ 'U nj nj /-iH ■t/ 1—1 

i — l — i ^ 2 , 2+1 ^ ^2+1,2] 5 

(Bl) 

as schematically shown in Fig. 0 Here = [E— < 
i\H\i >]“^ is the bare Green’s function of the isolated 
ith principal layer, with the superscript 0 indicating it 
is defined as if all other layers are deleted. The matrix 
T-Li^i-i =< i\'H\i — 1 > contains all the Hamiltonian ma¬ 
trix elements connecting sites in the layer *—1 to the layer 
i. Similarly Gf_ ^ =< i — 11 G^ | * — 1 > is a matrix defined 
on the {i — l)th principal layer, where G^ is the exact 
Green’s function of a subsystem of our original lattice 
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with all layers to the right of the {i — l)th layer deleted, 
as shown in Fig. The superscript “L” here means 
that this subsystem, including a left lead, is extended to 
the X = —oo. Since the (z — l)th layer is the surface of 
this subsystem, we will call Gf_i the left surface Green’s 
function. Similarly =< i + lIG-^lz + 1 > is another 
surface Green’s function of a subsystem of our original 
lattice with all the layers to the left of the {i + l)th layer 
deleted. Once Gf_j,Gf^^ are known, Gi can be com¬ 
puted immediately from Eq. |B1| 

The central task is then to compute Gfl^ and G^^. 
This can be done recursively. Take Gf_i as an exam¬ 
ple. We start with the leftmost layer i = 1. There our 
central system is connected to a semi-infinite lead, which 
contains infinite number of layers of the same width M, 
numbered by i = ..., —2, —1, 0. We denote this left lead’s 
surface Green’s function as G^, whose computation will 
be presented in the following appendix subsection |B 1| 
Then we add the i = 1st layer of our central system, 
but not other layers, to this lead so that we get a new 
semi-infinite stripe. This new stripe has a new surface 
Green’s function denoted as Gf, which can be computed 
from Gf' by 

Gi = [G?”' - ^1.0 G^ (B2) 

where Hifi connects sites in the surface layer z = 0 of the 
left lead to the z = 1st layer of our central system. Simi¬ 
larly we can repeat this process by adding one more layer 
of our central system to the semi-infinite stripe each time, 
and build up the whole system. In general at an inter¬ 
mediate stage, we may have a semi-infinite stripe, whose 
surface is, say, the (z — 2)th layer with a surface Green’s 
function Gf_ 2 - Then the (z — l)th layer is connected to 
that stripe to form a new semi-infinite system, which has 
a new surface Green’s function Gfl^. And Gf_i can be 
calculated from Gf _2 by the following recursive relation 


Similarly the right surface Green’s function Gf^^ can 
be computed from G| 5|_2 via 

Gf+i = [G°+r' - H.+i,.+ 2 Gf+2H,+2,,+i]-i (B4) 

This recursive relation starts with G £^2 the rightmost 
layer z = L/2 of our central system, where it is connected 
to another semi-infinite stripe lead extended to x = oo. 
Note that the central system has only L/2 principal layers 
because each layer contains two columns of the sites, and 
there are only L columns in total. The layers in this right 
lead are numbered by z = L/2-|- 1, L/2 + 2,.... We denote 
the right lead’s surface Green’s function as Gf. Then 
G ^/2 can be computed from G^ by 

Gl/2 = [Gl/ 2 ~ '^L/2,L/2+l Gf 'Hl/2+1,L/2] ^ (B5) 

where 'Hl/ 2 ,l/ 2 +i connects our central system to the 
right lead and contains t,t',t" only. 

1. Surface Green’s function ,G^ of the leads 

Gf',Gf can be computed by solving a self-consistent 
2x2 matrix equation. We now give a detail discussion on 
how to compute Gf', but only briefly mention the final 
results for Gf at the end. 

The left lead Hamiltonian contains only the hopping 
parameters 

0 M 

Hlend = —t[c|_i jCjj- -|- c| 

i— — OCl j—1 

~ t [cl-2,jG,j + A ~ LcJ jCij} 

(B6) 


Gf_i = [G°_i " -H^-l,^-2Gt2n^-2,^-l]-^ (B3) 

This is schematically illustrated in Fig. 


To be compatible with our central system Hamiltonian, 
which contains superconductivity, our lead Hamiltonian 
should have both an electron part and a hole part so that 
the full Hamiltonian Hiead is 





FIG. 8. Schematic diagram for the Gf_i computation. The 
sites enclosed by the blue dashed rectangle belong to the (z — 
l)th principal layer, which is also the surface layer of this 
semi-infinite stripe. 


^lead 


-fflead 0 \ 

0 —H\ea,d J 


(B7) 


Correspondingly the surface Green’s function takes a 
block diagonal form 


g,(l;+) 


[E+ - Lllead]-^ 

0 


0 ^ 

[E+ + Lliead]-^ J 


(B8) 


where for brevity we have introduced L+ = E + i6. We 
will denote the two diagonal terms as Gee = [E~^ — 
Hiead]~^ and Gth = [-G+ -b i7iead]“^- Apparently 
Ghh can be obtained from Gee by simple substitutions: 

^ Therefore we only 

need to discuss how to compute Gee- 
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Because of the periodic boundary condition along the 
y—direction, we can decompose Gee{E'^) into different 
momentum ky channels 


Gee{E+) ><XkJ g{ky,E+) (B9) 


with \xky >= ^7^ IJ > . and fcy = ^ with n = 

1, 2, 3, M. Each channel is described by a semi-infinite 
one dimensional chain effective Hamiltonian Heg^ky) =< 


XkjHieadlXky >, given by 
0 

Hee{ky) = ^ {{—21 cosky — 21" cos2ky — ^) clCi 

i— — oo 

+ [(-t -I- 2t' cosky) c-Ci_i - t"c-Ci _2 + h.c.]}. (BIO) 

And g{ky, E+) is the corresponding surface Green’s func¬ 
tion of this one dimensional chain. 

To compute g{ky,E'^) we group every two adjacent 
cites {c\_^ , cj) of the one dimensional chain together into 
a cell, indexed by the cell number n, so that Hee{ky) can 
be rewritten in a form such that direct couplings exist 
only between two nearest neighboring cells 


0 

Hes{ky) = ( C2„_i 

n— — oo 

+ X! ( 4n-l 

n— — oo 



—t" —t + 2t' cos ky 

0 -t" 


+h-c. 

C2n-2 J 


^2n 


—2t cos ky — 2t” cos 2 ky — fj, —t + 2t' cos ky 

—t + 2t'cosky —2t cosky — 2t" cos2ky 




C2n-1 

C2n 


(Bll) 


Since g{ky,E^) is a surface Green’s function, it s hould 
satisfy the same recursive relation given in Eq. (B3), 
which is rewritten here as 


1-1 


g — [^0 ~ [^^eff]o.-l G_i [iJeff]-l,o] 


(B12) 


The only difference from there is now all the matrix ele¬ 


ments are defined between different cells instead of layers. 
Eor clarity we have suppressed the ky and E+ depen¬ 
dence of all the quantities in this equation. Gjj is the 
bare Green’s function of the isolated single cell n = 0. 
Because each cell contains two sites, Gq is a 2 x 2matrix, 
given by 


G 


1-1 


= E+- [Efeff]o,o = E+- 


—2t cos ky — 2t” cos 2ky — g 
—t + 2t' cos ky 


—t + 2t' cos ky 
—2t cos ky — 2t'' cos 2ky — /i 


(B13) 


Similarly the effective hopping matrices between the 
cell n = 0 and cell n = — 1 can be read off directly from 

Eq. (IbTT) 


[^^eff]o,-l — 


—t" —t + 2t' cos ky 

0 -t" 


(B14) 

(B15) 


By definition G^i in Eq. (B12) is the surface Green’s 
function of the same chain but with the cell n = 0 deleted. 
However, since the chain is semi-infinite, deleting the sur¬ 
face cell only gives another identical semi-infinite chain. 
Therefore G^^ should be the same as g. Then Eq. (B12) 
becomes a self-consistent equation of g as 




£1+ -I- 2t cos ky + 2t'' cos 2ky + g 


t — 2t' cos ky 
—t” —t-\- 2t' cos ky 


0 


-t" 


t — 2t' cos ky 
E+ -I- 2t cos ky + 2t" cos 2ky -\- g 


-t" 0 

—t + 2t’ cos ky —t" 


(B16) 


With this 2x2 matrix equation, for each ky, we solve for g numerically by iterations until the results converge. Then 
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the computed g{ky, E) is substituted back into Eq. (B9) 
of G'ee(E+) to get . 

Similar derivations can be carried out for the right lead 
Green’s function Gf. It turns out Gf = (Gf')"'", where T 
is the transpose operation. This result is a manifestation 
of the fact that the two semi-infinite leads can be con¬ 
nected to each other by a reflection symmetry operation 
along the x—direction. 


Appendix C: The ansatz 

71 


The pairing amplitude on the bond, that connects two 
nearest neighboring sites and r^-, is calculated by the 
following ansatz: 


|A,,| = A 


TeS 


with Teff given by 




(Cl) 


eff 






(C2) 


where = 
r j+rj 
2 


I ri-l-rj 


center 


- is the distance from the bond 
to the nth vortex center R„, q is some pos¬ 
itive number, and Ny is the total number of vortices. 

If we consider a special case that there is only one vor¬ 
tex, for instance the nth vortex, then Eq. (C2) is reduced 
to Teff = r„, and 


In other words, we can define the pairing amplitude |A„| 
for the case when only the nth vortex is present as follows 

+ r 


so that |Aij| = |A„|. 

When more than the nth vortex is present, \Aij \ should 
become smaller than |A„|. This requires rgff < r„ be- 
causejAijl is an increasing function of r^g, as seen in 
Eq. (Cl I. We sum the contributions from each vortex to 
|Aij| simply by adding the gth inverse moment(g > 0) 
of al l r„ together to define an effective distance VeS as in 
Eq. (C2). Using the gth inverse moment, instead of the 
gth moment guarantees that r^g < r„ when there is more 
than one vortex. Furthermore it ensures that the terms 
(-5-)? with small r„ on the right hand side of Eq. (C2) 
contribute more significantly than those with larger r„. 
This is consistent with the physical intuition that vor¬ 
tices nearby are more important in determining rgff, and 
therefore |Aij|, than those that are far away. Also when 
there are more vortices present, Ny becomes larger and 
the resultant |Aij| from Eq. (C2) becomes smaller. This 
again agrees with our expectation. 


The I Ay I defined above increases monotonically with 
the parameter q for a given vortices configuration. To see 
this we only need to show r^g increases with q. For that 
purpose we can rewrite Eq. ( |C2[ ) as follows 

loglAik = liog{i+^'(!Aill).} (C5) 
ref[ q n 


where rmin = inin{r„} is the distance between the closest 
vortex and the bond, and the prime sign in the summa¬ 
tion means this closest vortex is excluded. The right 
hand side of Eq. (C5) is a monotonic decreasing function 
of q because in the summation each < 1. Therefore 
Teff increases monotonically with q, so does |Ay |. 



FIG. 9. Oscillation spectrum of the DOS for the two fold 
DDW order model with different q values but a fixed A. Note 
that the vertical axis scale for g = 1 is different from those for 
other q values. The two vertical dashed lines mark the two 
frequencies = 525T, Fy = 966T. 

An appropriate value of q can not be determined with¬ 
out solving the whole problem self-consistently, therefore 
we performed simulations for different q to see if our con¬ 
clusions depend on q or not. One example data of the 
oscillation spectrum for the two-fold DDW order case is 
shown in Fig.|^ From this figure, we observe that the os¬ 
cillation amplitude decreases as q is increased from q = 1 
to g = 3. This is consistent with the analyses that |Ay| 
is a monotonic increasing function of q, since larger q 
gives larger |Ay |, which means stronger vortex scatter¬ 
ing and therefore stronger suppression of the oscillation 
amplitudes. 

Although the oscillation amplitudes can depend signif¬ 
icantly on q, the oscillation frequencies remain unaffected 
by varying the q values, therefore the conclusion of On- 
sager’s relation being robust against the vortex scattering 
does not depend on the value of q. 


Appendix D: Check of the results in Ref.[7| 

We have checked the Fig.l and Fig.3 of Ref.[7], using 
the same parameter sets, and find our conclusions remain 
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the same. For both these two cases, the normal state, 
without magnetic field, can be described by the following 
Hamiltonian 

+^' X 

<*,j> «i,j» 

+ X * . (Dl) 

<i,j> i 


In this Hamiltonian the third term is a two-fold DDW 
order(or the staggered flux state order), and rjij = ±1 is 
again the local d—wave symmetry factor. 


1. First consider the Fig.l of Ref. [7]. We use the 
same parameters t = l,t' = 0.3t, Wq = 1.0t,/i = 
—0.949t. The normal state Fermi surface consists 
of four hole pockets with an area ( 2 ^/a)^ ~ 2.5% 
each. Fig. IT shows the computed DOS. We see 
there is no noticeable shift in the oscillation fre¬ 
quency when the vortex scattering is present. 



20 40 60 80 100 

1/B 


FIG. 10. DOS oscillation for t = l,t' = 0.3t, VFo = 1.0t,p = 
— 0.949t. The unit for the field B is ? with 'Fq = hc/e the 
full flux quantum and a the lattice spacing. And the DOS 
unit is states/t . In the legends the “normal” means A = 0. 
The lattice size is L = 2000, M = 80, and in the Eq. ( |C2[ ) of 
Teff, q = 1 rather than g = 2 has been chosen here. 


2. Then consider the Fig.3 of Ref. [7]. In this case, the 
normal state does not have DDW, so Wq = 0. For 
t = l,t' = 0.14t,/i = —2.267t, the obtained Fermi 
surface contains only a large hole pocket with an 
area ~ 14% at the Brillouin zone center. 

Fig, [ill shows the corresponding DOS results. We 
see the oscillation amplitude gets heavily damped 
as A increases. Moreover, a small frequency shift 
SF/F « 2% becomes noticeable. However, this is 
different from a large 30% shift found in Ref. [7]. 
Also this 2% shift does not contradict our previ¬ 
ous conclusion of no noticeable frequency shift. Be¬ 
cause the shift here is obtained at magnetic fields 
that are larger than the experimentally applied 


fields(^ SOT) by an order of magnitude. In Fig. 11 


g = 10 corresponds to B = « 450T, since 

$o/27ra2 « 4500T if we take a = 3.83A for YBCO. 



1/B 


FIG. 11. DOS oscillation for t = l,t' = 0.14t, IFo = 0,/r = 
— 2.267t. In this simulation the system size is L = 1000, M = 
100 . 
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